package top.saodisheng;

/**
 * Description:
 * 求解矛盾方程最小二乘解应用模块
 * @author 扫地生_saodisheng
 */
public class CLEG_LSM {
    static double[][] PC;      // 矛盾方程组的系数矩阵
    static double[][] PA;      // 正规方程的增广矩阵
    static double[] PY;            // 矛盾方程组的常数向量
    static int PL;      // 系数方程的行数
    static int PM;      // 系数方程的列数
    static int PN;      // 正规方程增广矩阵的列数

    // 生成正规方程的增广矩阵
    static void FormNormalEquation() {
        double temp;
        // 计算系数矩阵与系数矩阵转置矩阵的乘积
        for (int i = 0; i < PM; i++) {
            for (int j = i; j < PM; j++) {
                temp = 0;
                for (int k = 0; k < PL; k++) {
                    temp += PC[k][i] * PC[k][j];
                }
                PA[i][j] = temp;
                PA[j][i] = temp;
            }
        }
        // 计算系数矩阵PC与常数向量的PY的乘积
        for (int i = 0; i < PM; i++) {
            temp = 0;
            for (int k = 0; k < PL; k++) {
                temp += PC[k][i] * PY[k];
            }
            PA[i][PM] = temp;
        }
        return;
    }

    static int Jordan() {
        double x;
        for (int k = 0; k < PM; k++) {
            x = PA[k][k];       // 对应行的对角线元素
            // 归一化操作
            for (int j = k; j < PN; j++) {
                PA[k][j] /= x;
            }

            // 消元
            for(int i = 0; i < PM; i++){
                if(i == k) continue;
                x = PA[i][k];
                for(int j = k; j < PN; j++)
                    PA[i][j] -= x * PA[k][j];
            }
        }
        return 0;
    }

    // 传入正规方程的解
    static double Objective(double[] A){
        double OBJ = 0.0;
        double temp;    // 保存孩子的预测身高
        for(int i = 0; i < PL; i++){
            temp = 0.0;
            for(int j = 0; j < PM; j++)
                temp += A[j] * PC[i][j];
            temp -= PY[i];      // 预测身高和真实身高的差值
            OBJ += temp * temp;     // 样本误差平方和
        }
        return OBJ;
    }


    /**
     * 把主程序传过来的参数复制到本模块的公共变量中
     * @param C ：矛盾方程组的系数矩阵
     * @param Y ：矛盾方程组的常数向量
     * @param A ： 问题的解向量
     * @param L ：系数矩阵的行数
     * @param M ：系数矩阵的列数
     * @return  ：返回结果码，0表示正确，-1表示样本数量过多或待求解的未知数过多
     */
    static double SolutionCLEG(double[][] C, double[] Y, double[] A, int L, int M) {
        int ErrCode = 0;
        if (M > 50) return -1;
        PL = L;     // 系数矩阵的行数
        PM = M;     // 系数矩阵的列数
        PN = M + 1; // 正规方程增广矩阵的列数
        // 初始化
        PC = new double[L][M];
        PA = new double[L][PN];
        PY = new double[L];

        // 将系数矩阵复制到本模块的公共变量PC中
        // 将常数向量复制到本地的功能模块公共变量PY中
        for (int i = 0; i < L; i++) {
            for (int j = 0; j < M; j++) {
                PC[i][j] = C[i][j];
            }
            PY[i] = Y[i];
        }

        // 生成正规方程的增广矩阵
        FormNormalEquation();
        // 调用Jordan()函数求解正规方程
        ErrCode = Jordan();
        // 无解
        if (ErrCode > 0) return -3;
        // 把正规方程的解放入数组A中
        for (int k = 0; k < PM; k++) {
            A[k] = PA[k][PN - 1];
        }

        // 利用Objective()函数计算目标函数的误差
        return Objective(A);
    }

}
